Replica Symmetry Breaking Condition Exposed 
by Random Matrix Calculation of Landscape Complexity 



Yan V. Fyodorov and Ian Williams 



School of Mathematical Sciences, University of Nottingham, Nottingham NG72RD, England 



We start with a rather detailed, general discussion of recent results of the replica 
approach to statistical mechanics of a single classical particle placed in a random 
N(':$> l)-dimensional Gaussian landscape and confined by a spherically symmetric 
potential suitably growing at infinity. Then we employ random matrix methods 
to calculate the density of stationary points, as well as minima, of the associated 
energy surface. This is used to show that for a generic smooth, concave confining 
potentials the condition of the zero-temperature replica symmetry breaking coin- 
cides with one signalling that both mean total number of stationary points in the 
energy landscape, and the mean number of minima are exponential in N. For such 
systems the (annealed) complexity of minima vanishes cubically when approach- 
ing the critical confinement, whereas the cumulative annealed complexity vanishes 
quadratically. Different behaviour reported in our earlier short communication [Y.V. 
Fyodorov et al., JETP Lett. 85, 261, (2007)] was due to non-analyticity of the hard- 
wall confinement potential. Finally, for the simplest case of parabolic confinement 
we investigate how the complexity depends on the index of stationary points. In 
particular, we show that in the vicinity of critical confinement the saddle-points 
with a positive annealed complexity must be close to minima, as they must have a 
vanishing fraction of negative eigenvalues in the Hessian. 

1 Introduction and General Discussion 

After more than three decades of intensive research, understanding statistical mechan- 
ics of disordered systems still remains a considerable challenge to both theoretical and 
mathematical physics community. In the situations where an interplay between thermal 
fluctuations and those due to quenched disorder is essential, one of the central problems 
is to find the averaged value of the free energy of the system as a function of temperature 
T. Namely, given an energy function which assigns a random value '^^(x) to every point 
X G S of a given configuration space S, the task is to calculate 



Here and henceforth the angular brackets stand for the disorder averaging, (3 = 1/T is 
the inverse temperature, and Zp is known as the partition function. In the case when 
the configuration space consists of discrete points (e.g. the hypercube S = { — 1, 1}^) the 
integration is replaced by an appropriate summation. 



Abstract 



F 





1 



Confronted with the problem of performing such a calculation explicitly, theoretical 
physicists almost invariably choose to employ the so-called replica trick - a powerful, yet 
heuristic way of extracting the above average from the positive integer moments of the 
partition function. The trick amounts to exploitation of the formal identity 



Indeed, assuming the distribution of random energies 7i(x) to be Gaussian with the pre- 
scribed covariance in the configuration space, the moments can be readily calculated for 
all values n = 1, 2, .... Unfortunately, for all systems of interest those moments grow too 
fast with n to allow a unique restoration of the probability distribution of the partition 
function, hence prohibiting the honest evaluation of the expectation value of the loga- 
rithm. Yet, physicists managed to employ their fine intuition and developed a refined 
and efficient heuristic machinery of performing a mathematically ill-defined continuation 
n — i> 0. In this way they were able to reconstruct highly nontrivial patterns of behaviour 
typical for disordered systems, the most intricate scenario being based on the notion of 
Replica Symmetry Breaking (RSB). According to that scenario first developed by Parisi 
in the context of the so-called mean-field models of spin glasses, see [Ij as the standard 
reference and p] for a recent account, below a certain transition temperature (frequently 
referred to as the de-Almeida-Thouless [3] temperature) the Boltzmann-Gibbs measure 
on the configuration space is effectively decomposed into a multitude of the so-called 
"equilibrium states" (this phenomenon is frequently referred to as an ergodicity break- 
ing). The number of those states is argued to grow exponentially with the number of 
effective degrees of freedom, and they are believed to be organized in the configuration 
space in an intricate hierarchical way. Essentially, the structure is thought to represent a 
picture of valleys within valleys within valleys, etc., in an effective random (free) energy 
landscape [1] arising in such models. As that unusual picture arouse via the application 
of a somewhat problematic replica trick, it is important to mention that a few rather 
nontrivial results emerging in the framework of the Parisi scenario were recently recov- 
ered by rigorous mathematical methods. These important developments are mainly due 
to recent seminal works by Talagrand E], based on earlier results by Guerra [7], see 
also interesting works by Aizenman, Sims and Starr [S]. In particular, Talagrand was 
able to demonstrate that the equilibrium free energy emerging naturally in the Parisi 
scheme of RSB is indeed the correct thermodynamic limit of the free energy, both for the 
paradigmatic Sherrington-Kirkpatrick model [6], and for the so-called spherical model of 
spin glasses 

One of the simplest, yet nontrivial representatives of a disordered system is arguably 
a single classical particle placed in a superposition of a random Gaussian V{x.) and a 
deterministic confining potential V^on(x), with x G S = M^. It turned out to be a 
surprisingly rich model, characterised by a non-trivial dynamical behaviour as well as 
interesting thermodynamics. Works by Mezard and Parisi [9], and Engel [TU] used the 
replica trick to calculate the free energy of such a system with the simplest, parabolic 
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confinement V^o„(x) = ^fix."^ , fj, > and after furtlier employing tlie so-called Gaussian 
Variational Ansatz (GVA) revealed the existence of a low-temperature phase with broken 
ergodicity. They were followed by Franz and Mezard [TTj and Cugliandolo and Le Doussal 
[T2] papers on the corresponding dynamics revealing long-time relaxation, aging, and other 
effects typical for glassy type of behaviour at low enough temperatures. The nature of 
the low-temperature phase was found to be very essentially dependent on the type of 
correlations in the random potential, specified via the covariance function chosen in the 
form ensuring stationarity and well-defined large- iV limit as 



with brackets standing for the averaging over the Gaussian potential distribution. Namely, 
if the covariance f{x) decayed to zero at large arguments, the description of the low tem- 
perature phase was found to require only the so-called one-step replica symmetry breaking 
(IRSB) Parisi pattern. In contrast, for the case of long-ranged correlated potentials with 
f{x) growing with x as a power- the full infinite-hierarchy Parisi scheme of replica 
symmetry breaking (FRSB) had to be used instead. 

Based on formal analogies with the Hartree-Fock method Mezard and Parisi [S] argued 
that GVA-based calculations should become exact in the limit of infinite spatial dimension 
N. In a recent paper |13] the replicated problem was reconsidered in much detail by an 
alternative method which directly exposed the degrees of freedom relevant in the limit 
N oo, and employed the Laplace (a.k.a. saddle-point) evaluation of the integrals. The 
method allowed also to perform calculations for any fixed sample size |xp < L = R\/N. 
The effective radius R < oo can be used as an additional control parameter, and the 
results in the limit R ^ oo and fixed > indeed reproduced those obtained by GVA 
in [9l Ho]. In addition, the research in [13] and the subsequent work [14j revealed the 
existence of a non-trivial class of random potentials overlooked in earlier papers[9],|T0], the 
simplest representative being the case of logarithmically growing correlations f{x) = /o — 
(7^ In (x + a^), with /o, g, a being some positive constants. Indeed, for such a choice (and its 
generalizations [T?]) the resulting free energy in the thermodynamic limit R ^ oo turns out 
to coincide precisely with that appearing in the celebrated Generalized Random Energy 
Model (GREM) by Derrida[T5l [16], which is, in particular, known to describe directed 
polymers on a tree with disordered potentials [T7]. Moreover, by comparing the results 
of [H] with earlier renormalization-group analysis of same system for finite dimensions 1 < 

< oo, see Carpentier and Le Doussal [T8], it was conjectured that essentially the same 
type of behavior should survive in finite-dimensional systems. Note, that understanding 
the behaviour of the model for finite N in full generality remains a challenging problem, 
see e.g. [2D] for a recent discussion. 

According to a general wisdom, many peculiar features of systems with quenched 
disorder as compared to their more ordered counterparts are largely due to a multitude 

-^To be more precise, at large separations (^[V (xi) — V (x2)]^^ cx (xi — X2)'^'^ , < 7 < 1. However 
one can easily satisfy oneself that the difference is immaterial for the free-energy calculations. 
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of nearly degenerate (metastable) states making the associated (free) energy landscape 
extremely corrugated, see e.g an insightful paper [19j trying to understand the gross struc- 
ture of the energy landscapes describing elastic manifolds in presence of a random pinning 
potential. In particular, for the paradigmatic mean-field models of spin glasses a consid- 
erable effort has gone into trying to understand what changes the intricate scenario of the 
replica symmetry breaking interpreted as ergodicity breaking (see our earlier discussion 
of the Parisi scenario) implies in the statistics of the associated free energy landscape (the 
so-called Thouless- Anderson-Palmer (TAP)[4J variational functional). For an instructive 
review of the recent efforts in this direction a reader may wish to look through [2T], or to 
read directly representative original works on this topic |22j-|27j. 

For the model of interest for us in this paper - a single particle in a random Gaussian 
potential - the necessity (and difficulty) of adequately taking into account multiple minima 
in the energy landscape is well-known since the early work of Engel|28] who pointed out 
a non-perturbative character of the problem. In fact, similar non-perturbative effects are 
believed to invalidate the famous "dimensional reduction" phenomenon in the random- 
field Ising model, see [29| [30] and references therein. It is therefore important to develop 
methods and tools of statistical characterization of relevant features of A^— dimensional 
random energy landscapes. Ideally, one would like to obtain a detailed knowledge on 
typical positions and depths of various minima and their mutual correlations. In such a 
general formulation the problem looks extremely challenging mathematically, although in 
the dimension A^ = 1 considerable progress can be achieved fST]. 

Obviously, the first meaningful step in this direction consists of counting the expected 
number of stationary points of different types (minima, maxima and saddles of various 
indices). The simplest, yet already a non-trivial problem of this sort is to find the mean 
number (i^s) of all stationary points, irrespective of their index. Assuming that the 
landscape is described by a sufficiently smooth random function Tioi N real variables x = 
(xi, xat) the problem amounts to finding all solutions of the simultaneous stationarity 
conditions dkH = for all = 1, A^, with dk standing for the partial derivative 
The total number ^s{D) of its stationary points in any spatial domain D is then given 
by #s(-D) = /^Ps(x)(ix, with Ps(x) being the corresponding density of the stationary 
points. The ensemble-averaged value of such a density can be found according to the 
multidimensional analogue of the so-called Kac-Rice formula: 



where 6{x) stands for the Dirac 6— function. For the original Kac and Rice papers in 
A^ = 1 see [33| [3l]. the multidimensional generalization can be found e.g. in [35] and [36] , 
and [37] contains further extensions and a proof in general setting. Note importance of 
keeping the modulus in (jlj), as omitting it would yield instead of the density the object 
related to the Euler characteristics of the surface, see [38] and discussions in [37] (see also 
the related identity (fTOj) below). 
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Similarly, if one is interested in counting only minima, the corresponding mean density 
can be written as 

(p™(x)) = (det {dl,,n)e {d^.^n) n'^(^^^)^ ■ (s) 

Here and henceforth we use the notation 6{x) for the Heaviside step-function, i.e. 6{x) = 1 
for a; > 0, and zero otherwise. The corresponding matrix ^—factor in ([5]) selects only 
stationary points with positive definite Hessians, which are minima. 

In the context of the problematic of glassy systems calculation of those quantities 
for the simplest case = 1 were presented in [32]. Actually, low- dimensional cases 

= (1,2,3) were known much earlier, starting from the classical papers by Longuet- 
Higgings on specular light reflection from a rough sea surface [39], and the work on laser 
speckle patterns by Weinrib and Halperin[30]. The case A^ = 3 was also considered in 
early work [H] in the context of behaviour of a quantum particle in a random potential, 
and recently in [12]. A useful summary of related efforts in various areas of applied 
mathematics and probability theory can be found in a recent book by Adler and Taylor [3 7J. 

Random landscapes relevant in statistical physics are however mainly high-dimensional. 
In particular, general interest in landscapes behind glass-forming mixtures of various types 
led to a few numerical investigations of the structure of their critical points, see e.g. [i3] - 
[IB] . Independently, interest in counting extrema of some random high-dimensional sur- 
faces arose also in string theory where an effective string landscape is believed to possess 
a huge number of possible minima, each of which could describe a potential universe, 
see [17] for further discussion, references and actual calculations. At the same time, no 
analytical results on numbers of critical points were available for high- dimensional Gaus- 
sian landscapes until very recently, before it was realized that obtaining them requires 
exploitation of well-developed random- matrix based techniques ^48j. For the A^— dimen- 
sional Gaussian random surface with parabolic confinement, i.e. for the function 

Hpar(x) = |x2 + V(x), ^>0 (6) 

the total expected number of stationary points calculated from the Kac-Rice density (jlj) 
turned out to be given by [18] : 

(#.) = ^(|det(/x/;v + ^)|) , (7) 

where H stands for N x N Hessian matrix of second derivatives of the Gaussian part 
of the potential: Hij = dfjVix.) , = 1, A^}, and In for the identity matrix, and 

we took into account that the distribution of the Hessian is position-independent due to 
translational invariance of the covariance function ([3]), see Eq.® below . We see that 
the problem basically amounts to evaluating ensemble average of the absolute value of 
the characteristic polynomial of a certain real symmetric random matrix, whose entries 
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are centered Gaussian variables due to the nature of the underlying potential V(x). The 
form of the covariance function Eq.Q implies after a simple calculation the following 
covariances of the entries Hij : 

/"(O) 

{HilHjm) = ^ [SijSlm + SimSlj + SilSjm] ■ (8) 

where here and henceforth in the paper the number of dashes indicates the number of 
derivatives taken. This allows one to write down the density of the joint probability 
distribution (JPD) of the matrix H explicitly as [4S] 



V{H)dH cx di/exp 



N 



4/"(0) 



Tr(i/n-^^fTr^^' 



(9) 



where dH = ni<j<j<Ar dHij and the proportionality constant can be easily found from 
the normalisation condition. Although such a JPD is invariant with respect to rotations 
H — s> 0~^H0 by orthogonal matrices O G 0{N) it is apparently different from the 
standard one typical for the so-called Gaussian Orthogonal Ensemble (GOE) of random 
symmetric matrices (the standard introduction into the Random Matrix Theory is [49j). 
A peculiarity of this ensemble is particularly manifest via the following identity that holds 
for an arbitrary fixed matrix A [12] : 



(^det{A + H)^ = detA. (10) 

In particular, this identity gives another evidence for importance of keeping the modulus 
in formulae like Eq.([7]), where omitting the modulus would trivially yield unity due to (fTOj) . 
An elegant general verification of (ITOl) directly from ([8]) can be found in [42j. Alternatively, 
in Appendix A we demonstrate its validity for real symmetric matrices by a specific 
method exploiting the very meaning of the matrix H as the Hessian matrix, as the latter 
is related to topological properties of the underlying gradient field. 

In spite of all the mentioned peculiarities, actual evaluation of the ensemble average in 
Eq.([7]) can be reduced to the standard random matrix calculation of the GOE density and 
as such can be performed in a closed form for any in terms of Hermitian polynomials 
[15] . This is important, as the same method actually proved to be useful in other similar 
problems, as e.g. the problem of determining the distribution of the maximum of a 
Gaussian random field, see [50] . 

Large asymptotics of the mean number of stationary points can be extracted from 
known behavior of Hermite polynomials. The explicit calculation revealed [18] that the 
replica symmetry breaking (interpreted in the standard way as ergodicity breaking) in 
the zero-temperature limit of the R = oo version of the model with parabolic confine- 
ment is accompanied by the emergence of an exponentially large mean total number of 
stationary points in the energy landscape. As common wisdom was to expect that the 
low-temperature thermodynamics should be dominated by minima rather than by the 



6 



totality of stationary points the issue called for further investigation, in particular in a 
general R < oo case. In two recent short communications [STl [52] the calculations of [48j 
were independently and essentially generalized in a way providing access to the densities 
of minima [32], as well as to stationary points with an arbitrary index [JI], and at a given 
value of the potential, hence energy. 

The results in [52] were given as a function of both parameters fi and R. That analysis 
revealed that generally the domain of existence of the glassy phase with broken ergodicity 
(at zero temperature T) turns out to be always associated with the existence of exponen- 
tially many stationary points in the energy landscape, but not necessarily exponentially 
many minima. In an attempt to extend those considerations to finite temperatures the 
authors also constructed a simple variational functional providing an upper bound on the 
true free energy of the /i > 0, i? = oo version of the problem. Surprisingly, counting 
stationary points in that simple-minded approximation was already enough to capture 
such a nontrivial feature as the precise position of the de-Almeida- Thouless line in the 
whole (/X, T) plane. 

The aim of the present paper is to describe the random-matrix evaluation of landscape 
complexities in a rather general setting. The calculations are performed for a broad class 
of Gaussian landscapes V(x) with monotonically increasing, concave confining potentials, 
i.e. for the energy functions 



assuming U'{z) > 0, U{z)" > 0, V2; > 0, with the first (confining) term chosen in the 
rotationally-invariant scaling form ensuring a well-defined large N limit of the model. 
Solving this problem will also be used as an opportunity to provide a rather detailed 
exposition of the techniques sketched in our earlier short communication |52j . We will be 
able to treat in considerable generality the smooth confinement case, and compare it with 
non-analytic (hard-wall) confinement studied in [52] (and independently by A Bray and 
D Dean in [5l]). 

Our main results can be briefly summarized as follows. Define the function 



We find that for a generic, smooth concave confining potentials U{z) with a continuous 
derivative fi{z) the condition ensuring that the mean total number of stationary points 
in the energy landscape is exponentially large in coincides with one ensuring that the 
same is true for the mean number of minima. Different behaviour of two complexities 
for a particle confined in a finite box reported in [STl [52] is shown to be a consequence 
of non-analyticity of the hard-wall confinement potential. For a smooth confinement the 
common condition of positivity of the corresponding landscape complexities (defined as 
the logarithms of the mean total number per degree of freedom A^) reads 
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As shown in Appendix D, such an inequahty is precisely the condition of instabihty of 
the zero-temperature rephca symmetric solution of the associated problem in statistical 
mechanic^, see Eg. ( 11521) . The corresponding annealed complexity of minima vanishes 
cubically when approaching the critical confinement (M?]) . whereas the cumulative an- 
nealed complexity vanishes quadratically with the distance to criticality. Finally, in the 
appendix E we further investigate, using the method due to Bray and Dean[51j the an- 
nealed complexities of stationary points at a fixed value of the index I (the number of 
the negative eigenvalues of the Hessian). The calculation is performed for the simplest 
case of parabolic confinement Eq.([n]), when fi{z) = const = /i, hence Eg. (1601) amounts to 
/i < /icr = \/7''(0)- The results reveal that the only stationary points with non-vanishing 
annealed complexity arising precisely at the critical value n = ficr are those for which the 
number of negative directions is not extensive: limTv^oo ^ = 0- we move inside the 
glassy phase away from the critical value = ficr, stationary points with an increasing 
range of indices start to have a positive complexity. 

In the Appendices A,B, and C we provide a few useful technical details, some of 
them are not immediate to find in the available literature. We comment on the geomet- 
ric features of the landscapes behind the curious formula ffTOl) in the Appendix A, and 
discuss the distribution of the diagonal element of the resolvent of a random GOE ma- 
trix in the Appendix B. In the Appendix C we provide a short overview of a powerful 
heuristic random- matrix technique due to Dean and Majumdar [56] based on the func- 
tional integration. We found this method indispensable when calculating the complexity 
of minima. Finally, a statistical-mechanics calculation used for extracting the point of 
zero-temperature replica symmetry breaking for the present class of models within the 
framework of the replica method is sketched in the Appendix D, and the calculation of 
complexity of stationary points with a given (extensive) index is sketched in Appendix E. 
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"Random Matrices and Polynomials: a tool to understand complexity". The authors are 
grateful to H.-J. Sommers for collaboration at early stages of the project (see papers |13], 
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^This result, in particular, gives a kind of aposteriori justification of the use of the mean value of 
stationary points/extrema as sensible characteristics of high-dimensional landscapes. Indeed, in the ideal 
world we would like to know the typical rather than the mean values, i.e. to replace "annealed" averages 
featuring in ([4]) and (O with the corresponding "quenched" ones. The latter would require performing 
the averaging of the logarithm of the number of stationary points. Unfortunately, the present level of art 
in the field makes such a quenched calculation hardly feasible. 
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2 Complexity of stationary points versus complexity 
of minima 



2.1 The density of stationary points for spherically-symmetric 
potentials 



Our goal is to study the density of stationary points of the function (ITTI) around position 
X, subject to the condition that the random part V{x.) takes the prescribed value V. This 
is given by a variant of the Kac-Rice formula (jlj): 



(p.(y,x, [U])) 



det 



6 



5[V - \/(x) 



(14) 



where we used self-explanatory short-hand notations. 

For the choice of the confining potential in Eq. ( fTTl) we obviously have for the gradient 
vector and for the Hessian matrix 



dn 



dx. 



N 



dV 



dx. 



/i 



In + M(x) + 



N 



dxdx 



(15) 
(16) 



where we used the function fi{z) defined according to fll2p . as well a.s N x N rank-one 
matrix M(x) with entries expressed via components Xj, z = 1, . . . , of x as 



[M(x)] 



u" 



X 

2iV 



XiXj f 



N 



X 

iV 



(17) 



Note that for the parabolic confinement case Eq.([6]) obviously fi{z) = fi. We assume in 
the subsequent analysis that fi{z) is non-decreasing: f^'{z) > ,Wz > 0. 

We find it also convenient to collect all disorder-dependent factors in the function 



J^s{V,x,K, [U]) = {5[V-V{x)]S 



X /i 



+ 



dV 
dx 



6 



K 



dxdx 



and use the above expressions to rewrite the density Eq. ffT4l) in the form 



{ps{V,x, [U])) 



dK 



det 



X 



+ M(x) 



W,x,k,[U]). (19) 



In the rest of the calculations we assume that the covariance function f{x) in Eq.([3]) 
has finite values of the first two derivatives at the origin: |/'(0)| < oo and |/"(0)| < oo. 
Actually, from the point of view of statistical mechanics such a condition ensures the 
existence of a nontrivial phase transition at zero temperature, see Eq.f llSp . Note that 
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in the case of long-ranged potentials behaving like f{x) ~ x'^ with < 7 < 1 at large 
X this condition requires imposing some regularization at small arguments, e.g. f{x) = 
(a + xy, a > 0. 

Let us introduce the following notations for the subsequent use, remembering /'(O) < 

0: 



/^c. = /m, Rcr = v\f'm/no), 9' = no) 



/(O) 



> 0. 



(20) 



Our analysis of Eq. (fT8|) starts with introducing the Fourier integral representation for 
each of the delta-functional measures, - of scalar, vector, and matrix argument, corre- 
spondingly. This step facilitates performing the ensemble average explicitly as averaging 
exponentials containing terms linear with respect to Gaussian variables requires only the 
knowledge of their covariances. In the course of this procedure one has to exploit a few 
identities which follow from Eq.Q after simple calculations, as well as from the expres- 
sions Eq.([8]). Namely, 



{V\^)) = iV/(0), 



(9x (9x 



-/'(O), 



V 



dV 



dV d'^V \ 



(9x / \ (9x dyidx. / 



0, 



(21) 



Tr A 



N 



no) 



2TtA'+ (TtAY 



(22) 



(23) 



where a and A are an arbitrary A^— component vector and an x real symmetric 
matrix, respectively. 

After these steps one can easily integrate out the scalar and vector Fourier variables as 
they appear as simple quadratic terms in the exponential (the integration over the matrix 
Fourier variable requires a little bit more care, and is done at the next step). Thus, we 
arrive at the following expression: 



J^(V,x,K, [U]) oc exp 



2/'(0) 



2Nf{0) 



i{k) 



(24) 



where 1{K) stands for the integral over the remaining Fourier variables, that is over a 
real symmetric N x N matrix P, and is given by 



I(K) 



VI 



N 



^(TrP) \dP. (25) 



Here, and henceforth, we will systematically disregard various multiplicative constant 
factors for the sake of brevity. They can always be restored from the normalisation 
conditions whenever necessary. 
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To deal with the remaining integral over the matrix P in fl25|) in the most economic 
way it is convenient to employ first a Gaussian integral over an auxiliary scalar variable 
t, and in this way to "linearize" the last term in the exponential in Eg. (125 p : 



e 2iv 



N 
2^ 



(26) 



Substituting this identity back into fl25l) converts the matrix integral into a standard 
Gaussian one which can be immediately performed, yielding 



J^{V,x,K, [U]) oc exp 



2/'(o) r vA^yj 



X 



2iV/(0)j 



(27) 



X 



_ JV j2 N 

2 3^1 



-Tr[i'-(,i+^^)/. 



dt 



Before substituting this expression back into Eq. lUQl) . we find it convenient to rescale 
coordinates and the value of the potential as x — > y/Nx. and V — > NV, and to introduce 
two matrices H and Ai defined as 



and 



H = K-{gt+^V]I^ 



M = sIn + 2 x ® x^ , 



where we introduced a short-hand notation 

1 



(28) 
(29) 

(30) 



and the diadic product x ® X"^ stands for the (rank-one) N x N matrix with components 



As a result, we find that the mean density of stationary points Eq.f llOp is expressed in 
terms the rescaled variables as 

p(y,x, [U]) =KewN (x^) - ^) / dte-'^^'/'V^it^V, [U]) (31) 



v,{t,v,[u]) = ( det(i^ + ;w) 



GOE 



(32) 



where A/^ is the required normalisation constant and {■ ■ ■)goe denotes the average with 
respect to the standard GOE measure V{H)dH oc exp(— TtH'^/4) dH. Now we can 
use the rank-one character of the second term in fl2^ which implies the following identity 



det{H + M) 



1 + 2 



det if + si 



(33) 
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where 

GhM = . ^ . X (34) 
H + sI 

is a diagonal element of the resolvent of the random matrix H + s I. Note that with the 
parameter s being real as required by the identity (133!) . the resolvent can take arbitrary 
large values when an eigenvalue of the random matrix H occurs close enough to —s. As 
in the large N limit the eigenvalues of GOE matrices fill in densely the interval (—2, 2)[49J 
we can expect the resolvent to retain strong fluctuations for \s\ < 2. Outside that interval 
however we may expect self-averaging of the resolvent. The explicit calculation performed 
in the Appendix B confirms this intuition, and establishes the magnitude of the fluctua- 
tions inside the spectrum. Namely, we show that for N oo the distribution function 
of GHi'x) tends to the Cauchy law centered around the mean value {Gh{'x))goe = l-^x^ 
and with the widths T = (|G'//(x) - {Gh{^))goe\)goe = - for |s| < 2, and T = 
otherwise. We therefore conclude that the first factor in Eg. (1331) is typically of the order of 
unity at every realisation of the disorder. Recalling that only factors of order exp{0(A^)} 
in the mean density of stationary points are relevant for calculating the complexity, we can 
safely disregard the mentioned factor in the rest of the calculation by replacing Eg. (1321) 
with 

V,it,V,[U]) = (detiH + si) ) (35) 

and remembering the relation ( l30l) of s to parameters of the problem. 

For the asymptotic analysis of complexity we further use the result proved in ^8] . 
It essentially claims that the expectation of the determinant in the right-hand side of 
Eg. (!35|) can be represented in the large N limit as 

V,{t,V,[U]) (xe^vN^{s) (36) 

where the function $(s) = $(— s) is given explicitly by 

.2 v.^/^rzri , fs-^/s^ 



+ In 



(37) 



The method used in [H] was based on relating Vgit, V, [U]) to the mean eigenvalue density 
of GOE, and subseguent lengthy saddle-point analysis of an integral representation for 
that density. Perhaps, a shorter way to understand the above asymptotics is to notice 
that actually, 

^s)= / \n\s + \\ psciX) d\, p,,{\) = —Va^ (38) 

J-2 2tx 

where the integral is understood in the sense of the principal value. The above formula 
can be verified by applying ideas from statistical mechanics to the evaluation of -Ds(t, V^). 
A mathematically rigorous method of analysing this type of problems can be found in 
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a detailed paper by Boutet de Monvel, Pastur and Shcherbina [55]. An alternative, and 
very transparent way of heuristic asymptotic analysis based on the concept of a functional 
integral was proposed in a recent insightful work by Dean and Majumdar (SB]. We provide 
an overview of the Dean-Majumdar approach in the Appendix C. 

We thus arrive to the following expression for the density of stationary points 

/oo 
■QO 

valid in the large- limit up to factors of order of unity. 



p(\/,x, [U])^Kexp{N 



2/'(0) 



2/(0) 



2.2 The cumulative complexity of stationary points for a smooth 
confinement 

The mean total number of the stationary points for a given shape of the spherically- 
symmetric confining potential is obtained by further integrating this density over the 
position x and over the values of the potential V: 

(#,) = / rfx r p{V,^,[U])dV (40) 

It is natural to use the fact that the integrand is spherically symmetric by passing to 
the integration over the radial variable q = and the angular coordinates fi, so that 
fix oc q^^~'^^/'^dqdVt. As a result, we have 

(#^) oc r e^^(^) ^^(g), ^(g) = H H e'^'^t^''^'*] (41) 

Jo Q J -oo J -oo 



5(^,.,t)^ ^4 -$(.), £(,)^^A.) + ^ln., (42) 



(43) 



where 



with (cf. ([30])) 

^ _ 1 

and $(s) given by Eq. fl371) . As usual, we disregarded in fHTj) all constant proportionality 
factors which will be restored at a later stage from the normalisation condition. 

For our analysis of this expression in the limit ^ 1 we find it convenient to evaluate 
first the integrals over t, V with the Laplace method. We find T{q) oc exp [— A^iS(K, t*)], 
with K and being the values minimizing iS(V,g, t) for a given g, and satisfying the 
system of equations: 

_ j_rf$(£), 
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where according to Eq. (H3|) 



s* = —Hq) + 9t. + ^ K] = — [/i(g) + ^ K] , (45) 

/icr /(O) Her / (0) 



the second equahty following from the obvious relation = gV^,/ /'(O) implied by Eq.(l44 
and the definition of in Eq.(l20l). Moreover, for s* > Eq.( l371) imphes that 

d$(s), 1 



ds 



s,-9{s,-2),/sl-4: (46) 



Assuming first < < 2, the first of the equations (jUj) together with Eqs . 045146 1) 
immediately yields 

K=/^(g)^^^, = 2/i(g)— , andt* = /i(g)^ (47) 

l^cr l^cr ^^cr 

Obviously, this solution is compatible with our assumption < s,, < 2 only as long as 
< /i(q') < /icr- On the other hand, assuming > 2 the first of the equations (jH]) solved 
together with Eq. fl46l) yields the relation 

which is consistent with our assumption. Substituting for from Eq. (H5l) one immediately 
finds ^ 

V; = ^, U = -^ for /i(g) >/ie,. (49) 

In fact, V* is nothing else but the most probable value of the potential V{'x) at a stationary 
point of the energy surface 7i(x) situated at the distance |x| = yNR from the origin. 
We see that K experiences a a drastic change in its behaviour at the values of radial 
parameter q at which the curve /i(g) crosses the value ^cr = \/ /"(O) and zero. To 
simplify our analysis we assume in the following that /i(g) > and /i'(Q') > 0, Vg > 0, i.e. 
that yLi(g) is a non-decreasing function of its argument. Then the equation yu(g) = [i^r has 
either a single solution q = q^. when /i(0) < /icr < fJ'ig oo), or no solution at all (this 
happens if ficr < /^(O), or /i^ > fi{oo). 

The corresponding values of S{V^,t^:,q) can be found by substituting (H71) and (H9i) 
into Eq. fH2l) . In particular, the expression fHHl) in the regime (jUj) implies the identity 



— ^y si — 2] = SO that $(s*) = ^ + | — In and we obtain after some simple 
algebra: 



5(K,t*,g) = <| '^v;^^-^- (50) 
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Substituting J-'{q) oc exp [— A^iS(K, t,,,, g)] back to the integral fHTj) . we see that the 
total number of stationary points is given by 

(#.) = M r ef *(^) ^ , (51) 
Jo Q 

where A/l encapsulates all the necessary normalisation factors, and 




^^MW + M^ + ing, for /i(g) < /ic 



To find Afs it suffices to consider the special case where ^ oo when only the minimum 
at the origin survives and hence 1. This condition fixes with exponential in N 

accuracy 

AT, ~ e-^i'^^- (53) 

where the scale Rcr was defined in Eq.( l20l) . Taking this factor into account, and assuming 
that the equation fi{q) = ficr has a single solution g = g*, we split the integration range 
in flSTj) into two natural pieces < g < g^, and g > g*. Then the exponential in N part of 
the mean number of stationary points can be found from 



Jv,T, dp ]v,T,_ i„\ dq 



/ e^*<('^)^+ / eT*>('^)^, (54) 

where remembering f'{0) = —filr^cr defined 

^>(g) = 1 - w + \nw, where w = j^^, for fi{q) > ficr = /i(g*) 



^<(g) = (l - /!;) + In , for /i(g) < /i,, = /i(g,) 



(55) 



As the maximum of the function \l'>(g) = 1 — w + lnw is achieved at = 1 and equal 
to zero we have ^'>(g) < 0, Vg > g*, and hence only the first integral in flMl) can yield a 
positive complexity 

= hm 1 In (#,) = hm i- In f e^*<(^) ^] (56) 

To extract the large- asymptotics of the latter integral we calculate 
rf^<(g) dw fRl \ 1 / fi\q)\ 

— ( — +7 ( 1-^^ ' M?) < /i(g*) = /icr • (57) 



(ig dq \ q J q \ Z^; 
The analysis below proceeds separately for two cases: < g* < R^^ and < R^^ < g* 
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1. Let us first consider < R^^, with the goal to demonstrate that for such a choice 

dw 
dq 



there is no contribution to positive complexity. We have ^ > 0, so that the 



derivative in Eq- llFTI) is obviously positive: ^'^^^'^^ > 0, g < g*. Using that = 
w{q*) = < 1 we see \l/<(g) < \E'<(g*) = 1 — + In w,,, < and conclude that 
in such a case there is no exponentially large contribution to the integral, hence no 
positive complexity. 

2. Assume now < R^^ < g*. We shall demonstrate the existence of a point go G 
[Rlr, Q*] where the function \I^<(g) attains a positive maximum \l'<(go) > 0. This in 
turn will imply a positive complexity. 

To begin with, R^^ < g* implies fi{Rlr) < /^(o'*) = /^cr as the function /i(g) is non- 
decreasing. Then fl57j) implies ^^^^|g=R2^ = (^1 ~ ^^~^f^) ^ 0. Note that from 
fl55|l follows ^<(g < RIj.) < 0, and ^<(g = R^.) = 0. Therefore *<(g) > in some 
right vicinity of g = -R^r- 

On the other hand, we have already seen above that \I/<(g*) < and also ( 1571) 
implies 

c(^<(g) ^dw /g|, 
dq ''~'^* dq ^* \ g* 

Then by continuity there must exists a point g = gi G (-R^r' such that ^E'<(gi) = 0. 
Remembering \l/<(g = -R^^) = and positivity of ^E'<(g) to the immediate right of 
g = i?^^ we conclude on the existence of another point g = go G {R^^n 1i) where the 
function \&<(g) attains at least one positive maximum \l/<(go) > 0. 

Combining all these facts, we see that under these conditions Eq.f lS^ leads to the 
positive cumulative complexity of stationary points: 

> 0, (59) 



1 < . (5J 



^ 1 

= - max 



2 q€{Rl^,qt) /i^^ \ -^cr/ \^cr 

as long as fJ^iR^r) < fJ^cr- Taking into account definitions ( |20l) . the latter condition can be 
written in terms of the covariance function of the random potential as 

/^(-;^)<A/m. (60) 

This inequality is nothing else but precisely the condition of instability of the zero- 
temperature replica symmetric solution of the associated problem in statistical mechanics, 
see Eq.f ll52p in the Appendix D. 

For the simplest case of a parabolic confinement (jH]) with /i(g) = < fXcr we can easily 
find the cumulative complexity as an explicit function of fi. Indeed, the point where ^I/<(g) 

2 

attains its maximum is given by go = R^r^i 



- - 

2 



1 " 



1 1 ^ 
-1 - In — + — 

cr 



> 0, for yU < /icr (61) 
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in full agreement with the result reported earlier in [JH^. When approaching the critical 

value a = Ucr the cumulative complexity vanishes quadratically: T,s ^ (1 —) . 

Let us show that such type of critical behaviour is generic. Introducing the short-hand 
notations q = q/Rl^ and jl{q) = we rewrite \E'<(g) as 



f-l'cr 

;;\ r.2/~ 



^<(g) = (i-g)r(g) + ing. (62) 

The condition of criticality fl60|) is simply ^{E?^^) / = which amounts in new notations 
to = 1. Hence we introduce the parameter 5 = 1 — /i(l) which controls the distance 
to the criticality. Obviously, the interval i?^^ < q < q* shrinks to zero when 5 — 0, so 
everywhere in the critical region g — 1 = e ^ 1, and we can approximate the function 
by first two terms in Taylor expansion: jl{q) ^ +e/i'(l) = 1 — 5 + eji'{l). Substituting 
this to (!62|) we obtain 



vI/< (g = l + e) = (2<5 + 0(52))e-Q 



+ 2/i'(l) + 0(5) e2 + 0(e3). (63) 



To find the complexity, we should maximize this over e. To the leading order in 5 the 
maximum is attained at eo = ^/ (| + 2/i'(l)) and the corresponding complexity is given 
by 

1 (5^ 

= (g = 1 + eo) = rr , (64) 

s 2 < ^' 1 + 4/t'(l) ' ^ ' 

which indeed always vanishes quadratically at criticality. 

2.3 The density of minima for smooth spherically-symmetric 
potentials, and the associated complexity 

If one is interested in calculating the density of only minima, one has to perform the ma- 
nipulations identical to those in the first half of the preceding section leading to Eqs. fl31l) 
and (135!) . but using Eq.Q rather than (jl]) as a starting point. Repeating all the steps, 
one arrives then to 

p„(F,x, [f/])=Ar„expiv(^^/i2 (^')-^) / dte-''''''V^{t,V,[U])m 

with 

V^{t,V,[U])=D,r,{s) = Uet{H + si)e\det{H + si)\) , (66) 

where the variable s is given in terms of t and V by the same expression fl30p , and A/m is 
a relevant normalisation constant. As before we use the notation {. . .)goe for the average 
with respect to the standard GOE measure V{H) dH cx exp(— iV Trif^/4) dH. To extract 
large- behaviour of the function Dm{s) we employ again the Dean-Majumdar method 
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as explained in the Appendix C. In this way one finds that the function Dm{s) with the 
required accuracy is given asymptotically by 



exp|-^G^(s) + iVT^(s)} ifs<2, 
exp[Ar0(s)] if s > 2, 



Dm{s) OC 

where the calculation sketched in the Appendix C gives 



(67) 



Gmis) 



216 V / 6 



(6J 



so that Gm(2) = 0, and 0(s) is the same as expression (ITTI) for s > 2, i.e. 

o2 



sa/s2 - 4 



+ In 



(69) 



Explicit expression for the function Tm(s) can be easily found from the procedure described 
in the Appendix C, but is immaterial for our purposes, apart from the fact that continuity 
of V^s) at s = 2 and the property Gm(2) = implies Tm(2) = 0(2) = 1. 

Correspondingly, the density of minima in Eq. fl65|) is calculated as the sum of two 
contributions: 



P 



X 



exp 



^^ + NTUs)-—Gm{s)]e{2-s)dt, (70) 



P^^=^^expiv(^;.^(x^)-^) 



oo / jy^2 

exp — + N(P{s) ] 9{s -2)dt. (71) 



oo 



To shorten our analysis in the large- limit we recall that our final goal is rather to 
calculate the exponential in A^ contribution to the mean number of minima. As in this 
process we integrate the above density over the value of the potential V , and over the 
coordinates x in the sample, we can search simultaneously for the optimal values of 
variables t and V . Let us start with the asymptotic analysis of the second contribution, 
Eq. fITT]) . The shape of the integrand in this case coincides with that in the analysis of the 
density of stationary points for s > 2, hence the optimal values are given by (see (jUj)) 

K = ^^7^, = for /i(g) > yUcr • (72) 
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where as before we introduced g = x^. In particular, the most probable values of the 
potential V^(x) taken over all the stationary point, or over only minima coincide in this 
regime. The corresponding contribution to the density of minima is given by 



(2) ^ 

Pm exp — 



2/'(0) 



fi^ (q) - In 



(73) 



This is again the same as contribution of the density of all stationary points in the 
regime /x(g) > /icr- As we know from previous analysis , in this regime there can be no 
exponentially many stationary points, even less so minima. Similar consideration also 
makes it clear that the overall normalization factor Mm (common to both and p"m) 
in the final expression should be (within exponential accuracy) the same as one used for 
the total density of stationary points in Eg. (1551) . 

Now we consider the contribution coming from Eq. fl70p . which can be conveniently 
rewritten using the relation between t and s as 



Pm oc exp N 



2/'(0) 



p' (q) 



X 



exp-iV| — 




Tm{s) + —Gm( 



ds . (74) 



When N ^ oo, the integral over s will be obviously dominated by the value which 
minimizes the function Gm{s)- One can satisfy oneself that this minimum occurs precisely 
at the boundary of the integration region s = 2. Defining s = 2 — z we expand for small z 
as Gm{2 — z) ~ j^ + 0{z^), and after substituting this expansion back to Eq. flTij) integrate 
out z. This procedure yields an irrelevant pre-exponential factor, and using 0(2) = 1 the 
exponential in contribution to the density takes the form 



Pm <^ exp 



/'(o) 



pHq) 



/(O) 



2 - 



1 2^ 



(75) 



Next we find that this expression is maximized at the value of the potential given by 



m 

per 



2 - 



p{<i) 

per 



(76) 



and the value of the density pfn at this maximum is proportional to 



Pm^ oc exp 



N 



/'(O) 



p' (g) + 2 



2 - 



per 



(77) 



Now this expression can be used for extracting the exponential in contribution 
to the total mean number of minima, hence the corresponding complexity. This 
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amounts to multiplying the density (177|1 with the "volume" factor q^^~'^^^'^dq, integrating 
over the radial variable R in the range up to g = g^,, such that fi{q) < ficr for q < q^, and 
finally multiplying with the overall normalisation factor e~^^'^^='^. Taking the logarithm 
yields the complexity of minima 



lim — ln(#m)= lim — In 



- max "^miQ) 



dq 



where 



= 2 + In 



q fi^{q) 



2 - 



(78) 
(79) 

(80) 



Rl. fi, 

Let us first consider the simplest case of a parabolic confinement fi{q) = fi < ficr, 
see (|6]). The complexity of minima can be easily found as the function of ^/ ^cr- Again, 
the point go at which \I^m(g) attains its maximum is equal to the same value go = -R: 
which delivered earlier the maximum to the function \l/<;(i?), and 



'2 ^cr 



-3 - In — + 4 — 



for H < He 



It is evident that S^, ( — 



0, and it is easy to check that dT^m/dfi < 0, hence 



\ — — 1 ) > for any — < 1. Close to the critical value — = 1 the complexity 
~ i 

~ 3 



vanishes cubically: ~ | (^1 ~ ^ j • This is a faster decrease in comparison with the 
quadratic behaviour demonstrated by the cumulative complexity from (16T|) . In fact note 

that for any fi < ficr the inequality — = ^1 — > holds. 

In a general case /x'(g) > the point go at which the function \l/<(g) from (!55l) attains 
its maximum does not necessarily coincide with one which delivers the maximum to \I^m(g) 
from (IHOl) . Indeed, the two functions satisfy the following relation: 



^m(g) = ^<(g)-2 



1 - 



(82) 



which implies for derivatives 

d 



d 



^vE.^(g) = ^v,,(,)+4 



^3) 



Hence in the regime /i(g) < ficr interesting for us here ^\I'm(g) > at the point of 
maximum of the function ^<(g). 
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As discussed in the previous section, for q^, < R^^ the derivative ^\l/<(g) is strictly 
positive, and thus the derivative cannot be zero. This imphes that ^'^(g) can 

only have a maximum if the inequality > i?^^ holds, and the point of maximum must 
belong to the interval q G [R'^^,q^]. Below we shall prove that in fact the condition 
g* > RIj. is enough to ensure that complexity of minima is positive. Let us however 
note that the simple arguments used to prove the positivity of the maximum of \E's(i?), 
hence the positivity of the cumulative complexity Us, are not immediately applicable 
in the case of the complexity of minima. Indeed, although we have ^^l/m( 9)1^=^2^ = 

> 0, we still have (g = Rlr) < in view of the relation ( |82l) and 



~1? — 

r'cr 



\E'< (g = R^^) = 0. But we also have at the right end of the interval (q = g*) < and 
^^mlg* < 0, demonstrating impossibility to infer positivity of the maximum of \E'm('?) 
inside (i?^^,g,,) along these lines. 

To circumvent this difficulty, we introduce an auxiliary function \l'a(g, /i) of two real 
variables g and fl according to 

^,(g,/i) = 2 + lng-g/i2-[2-/i]2 . (84) 

Note that "^miQ) in flHOl) is obtained from the above function by replacing its arguments 
as g ^ g/-RL and u 

Considering values of "^aiQ, j^) we easily find that the function can take positive values 
only in a wedge-like region of the (g, jl) plane restricted by two boundary curves /i_ (g) 
and /i+(g), where 



^±(g) = ^4TI^' D{q)=A + {l + q){\n{q)-2) . (85) 

and one has to require -D(g) > for the existence of these curves. Since D{q = 1) = 
and ^^^^ = — 1 + Ing + i > 0, Vg > 0, we conclude that D{q > 1) > 0. This implies 
the region in between those two curves exists for g > 1, and for g — > 1 the two boundary 
curves approach each other and meet at the point fl±{q = 1) = 1 (see figured]). 

Now, the problem of finding the complexity of minima, Eq. fl78l) . for any monotonically 
increasing confinement function n[q) is equivalent to searching a maximum of the function 
\E'a(g, /i) along the corresponding curve /i(g) = n{q)/ in the plane (see the figure). The 
point q^, appearing in our analysis is nothing else but the point of intersection of that 
curve with the horizontal line /i = 1, and the condition g* > i?^^ just means that such 
an intersection happens at some value g* > 1. But every such curve /i(g) necessarily 
intersects the region between the boundary curves fi±{q), and therefore necessarily has 
a portion along which the values of \E'a(g, /i) > 0. This proves that the complexity of 
minima is always positive as long as the condition q^, > R"^^ holds. This is equivalent to 
n{Rlr) < fJ'cr which is precisely the replica symmetry breaking condition Eq. fl60|) . 

Let us finally demonstrate that for a generic smooth confining potential the complexity 
of minima vanishes always cubically at criticality, the type of critical behaviour we have 
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Figure 1: A representative function jl{q) which is monotonically increasing and satisfies 
the condition g,, > 1. These conditions are sufficient to ensure that the corresponding 
"^mici) has a positive maximum leading to a positive complexity of minima. 



already found in the simplest case of parabolic confinement. For this we follow the same 
steps as in the end of the preceding section. Using the parameter 5 = 1 — /2(1)^1 which 
controls the distance to the criticality, we approximate in the critical region g = 1 + e, e <^ 
1, and use for the function fl{q) its Taylor expansion. Truncating as before at first order: 
~ 1 — 5 + e/i'(l) and substituting this for ft in fl84|) . we obtain 



^„ (g = 1 + e) = -26^+ [26 + Afi'{l)6 - 6^)] e- 



i + 2^'(l) + 2(/i'(l))^ + 0(5) 



(86) 

To find the corresponding complexity, we should maximize this over e. To the leading 
order in 6 the maximum is attained at eo = 26/ (1 + 2/i'(l)). Substituting this value 
back to fl86|) we immediately find that all terms of the order of 6"^ cancel, and thus at 
criticality the complexity of minima behaves at least as = C 6"^. To prove that the 
coefficient C 7^ and to find its actual value in terms of requires going to 

the terms of order in the Taylor expansion for fl{q), and to the order in (186|) . This 
is a straightforward but boring exercise, and the resulting expression is not very elegant 
apart from the simplest case of the parabolic confinement when /i'(l) = = 0, and 

C = 1/3 as we have already found before. 

2.4 On anomalous critical behavior of complexities for hard-wall 
confining potentials. 

A fairly universal type of the critical behavior of both complexities and revealed in 
the preceding sections should be compared with the results reported recently in our short 
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communication [52] , and independently in [5T] . In both papers the landscape was confined 
in a finite box, chosen in [52] to be spherical of extent |x| < L = R\fN, < R < oo. On 
top of such a hard-wall confinement a harmonic confining potential was superimposed to 
allow a comparison with the results of statistical mechanics of the same problem studied 
by the replica trick in [13] . For the sake of clarity, we will consider below only pure hard- 
wall confinement, when the complexities for a given value of the radius R were found to 
be given by [51], 



(i?) = In ( ^ ) for R>Rcr, and {R) < otherwise , (87) 



and 



(i?) = — 1 + In ( -— ) for i? > Rm = eRcr, and (R) < otherwise . 



R, 



■cr 



with Rcr given by Eg. (1201) . Taking into account that the domain of zero-temperature 
replica symmetry breaking in this particular model is just given by i? > Rcrll3\, we 
see that apparently the type of behavior exemplified by Eqs. fl87j) . fl88|) is very different 
from what we have discussed earlier in this paper. Namely, although in this case the 
cumulative complexity is also positive everywhere in the phase with broken replica 
symmetry for R > R^r, the complexity of minima becomes positive only starting from 
a larger confining radius R^ = eRcr > Rcr- In other words, for the interval Rcr < R < Rm 
the broken ergodicity is not at all accompanied by the exponentially many minima in the 
energy landscape. Another peculiarity is that the complexity vanishes linearly rather 
than quadratically with the distance 6r = R/ Rc — I close to the ergodicity threshold. 

In the framework of the present approach the function fi{q) describing the hard- wall 
confinement with a radius R is formally f^{q) = for q < R^ and yu(g) = oo for q > R^. 
Such a form is apparently highly singular, and does not immediately fit into the analysis of 
the preceding sections. To circumvent this difficulty and to understand better the origin 
of the hard- wall behavior fl87|) . flHHj) within a more general framework we consider a family 
of non-smooth confining potentials of the form: 

r , for g < 

where > is a control parameter, and h{x) is an increasing, non-negative smooth 
concave function h{x) > 0,h'{x) > 0,h"{x) > for Va; > 1. We also assume h{l) = 
to ensure the continuity of fi{q), and thus existence of a single solution g = g^, of the 
equation /i(g) = /icr which as we know plays an important role in our analysis. We will be 
in particular interested in understanding the behavior of complexities T,s{R) and Sm(-R) 
in the limit /i ^ 1, where the confinement described by Eq. fl89|) should approach a hard- 
wall form. We also assume the confinement radius to satisfy R > Rcr, to allow for a 
positive complexity. 
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The cumulative complexity of stationary points for this class of potentials is given 
according to (159!) by si (i?) = |maXgg(i g^)\I'i'^^(_R, g) where 

^ ^ '^^ \\nq + fi^h{q/R^){l-q), for < q < q, ^ ' 

and we used as usual the scaled parameter /t = ^/ficr and the scaled confinement radius 
R = R/Rcr > 1 as well as g = q/Rlj., with the value of the parameter q^, being fixed by 
the condition jx^h{q^/ R^) = 1. Differentiating the expression (p7]) over q gives: 



^1 -^_~2 

dq^^-^'^- R2 ^ 



h'il){R'-l) , (91) 



which is clearly negative for large enough fi, except for the special case h{l) = h'{l) = 0. 

On the other hand, ^^lo<g<_R2 = > 0. We conclude that the function '$i^\R,q) for 

large enough /i has its local maximum precisely aX q = R^. Actually, this is the global 
maximum, as 

^U>R^ = -i {l + [2R' h'{q/R') + h"{q/R'){q - 1)] } < . (92) 

so that < ^|^^^2+ < for all R? < q < q^^. We immediately see that the cumulative 

complexity is then given simply by |lni?^. 

Let us now shortly discuss the complexity of minima along the same lines. The ana- 
logue of the expression (l90l) follows from ( IHOi) . and is given by 

_ ^ f lng-2, for < g < ^2 

*m - j i^^_2 + 4fi^h{q/R^)-fi\l + q)h{q/R^), for R^ < q < q, ^^^^ 

isily see that '^^'^ 



Differentiating the above expression, we again easily see that ^^|^^;j.2_,_ < for /i ^ 1, 



except for the case h{l) = h'{l) = 0. Since -^\g_^ji2_ is always positive, the function 

^m(-R; q) has its (global) maximum aX q = R^ and the complexity is indeed given by the 
expression (IHHj) . 

To investigate the behavior of both complexities in the remaining exceptional case 
h{l) = h'{l) = 0, we consider a particular example h{x) = [x — 1)^. This implies 

g* = -R^ ^1 + 4^ and for the case of cumulative complexity we are seeking to maximize 

^,(g) =lng + /i2(g7i?2_i)2(i_^) (94) 

over g in the interval i?^ < g < g*. For large ft it is actually more convenient to introduce 
a new variable y G [0, 1] via g = i?^(l + y/p) , and approximate \l/s(g) with 

^s{y) =\nie + y\\ - R^) + 0(1/^) , (95) 
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In view of the inequality > 1, the above expression attains its maximum at y = 0, 
and the cumulative complexity ^^(i?) is again given by the same expression (187|) as in all 
other cases of this family in the hard- wall limit /i ^ 1. 

Curiously enough, the complexity of minima in the case h{l) = h'{l) = 

shows a behavior different from flHHl) . Taking as before h{x) = {x — 1)^, we need this time 
to maximize 

^^(g) = Ing - 2 + Afl{q/R' - 1) - /^^(l + _ ^^2 ^gg^ 

over q in the interval < q < q* = (^1 + ■ Introducing again the variable y G [0, 1] 
via q = R?{1 + y / jl), and assuming yU ^ 1 we replace (1961) with 

^^{y) ^ lni?2 _ 2 + % + 0(1/^) -y\l + R') , (97) 

which has its maximum at ?/o = 2/(1 + R^) < 1. The corresponding complexity is given 
by 

4 

S„ = lni?2_2 + (98) 

which is indeed different from fl88p . In particular, the last term in fl98p ensures that the 
complexity of minima is positive everywhere in the phase with broken ergodicity R > 1, 
and vanishes linearly when approaching the critical value R> 1. 

Thus, we have demonstrated that for every curve in the family fl89|) the cumulative 
complexity in the limit fi oo is indeed given by S<,(i?) = ^InR^ = In R/R^r, in full 
agreement with the hard- wall confinement formula flHTl) . As to the complexity of minima, 
in the limit /i — cx3 it is generically given by (IHHI) . although the result may change if the 
hard-wall profile vanishes smooth enough when approaching the point of non-analyticity, 
q = R^. 



3 Conclusions and open questions 

Let us briefly summarized our findings. We have demonstrated that for a generic, smooth 
concave confining potentials with a continuous positive derivative the complexity ex- 
tracted from the mean number of totality of stationary points in the energy landscape is 
positive simultaneously with the complexity corresponding to the mean number of min- 
ima. The domain of parameters where those complexities are positive is precisely one 
where the zero-temperature limit of statistical mechanics in such a landscape requires for 
its description the concept of broken replica symmetry/broken ergodicity, Eq. fll52p . 

On the other hand, for a non-analytic (hard-wall) confinement the boundary of non- 
ergodic behaviour coincides in general with the domain of positive total complexity of sad- 
dle points but not necessarily the minima. Moreover, the above analysis clearly demon- 
strates that all peculiar features specific to the hard-wall confinement are due to the 
discontinuity of the derivative of the confining potential. 
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In the case of a smooth confinement the complexity of minima vanishes cubically when 
approaching the critical confinement (160!) . whereas the cumulative complexity vanishes 
quadratically with the distance to criticality. In the appendix E we follow the method by 
Bray and Dean [JT] and investigate the (annealed) complexities of stationary points at a 
fixed value of the number of the negative eigenvalues of the Hessian, the so-called index X of 
a critical point. Restricting the consideration to the simplest case of parabolic confinement 
Eq.Q we reveal that the only stationary points with non- vanishing complexity arising 
precisely at the critical confinement fi = ficr are those for which the number of negative 
directions stays of order of unity in the large— limit: a = limTv^oo ^ = 0- we 
move inside the glassy phase away from the critical value fi = ficr, stationary points with 
an increasing range of indices start to have a positive complexity. This is precisely the 
mechanism behind the change from cubic vanishing typical for the complexity of minima 
to the quadratic vanishing found for the total complexity. Whether the stationary points 
which appear precisely at the critical point are strictly minima, or may include admixture 
of saddle points with a few negative eigenvalues remains an interesting question deserving 
further, more elaborate, investigation. 

In our opinion, the results of present research suggest a few more immediate questions 
which are worth understanding. For example, the problem of investigating the magnitude 
of sample-to-sample fiuctuations of the number of stationary points is clearly very inter- 
esting, though technically challenging. Another natural extension would be to consider 
our model at finite temperatures. This project would require defining properly an ana- 
logue of the free energy (TAP-like, [4j) landscape for the present model, and recovering the 
whole transition line to the phase with broken ergodicity from counting the corresponding 
stationary points. Some steps in this direction were taken in [52] , but the present level of 
understanding remains far from being satisfactory. 

Finally, let us mention that the model considered in the present paper is very inti- 
mately related [13j to the much-studied spherical model of spin glasses, which essentially 
corresponds to a particular choice of the confining potential Vcon(x) = (5(x^ — N) which 
forces the vector x to span the sphere of radius V^. It is natural to expect that this 
rather singular limit can be approached in the present method by e.g. considering the 
confinement of the form U{z) = fi{z — 1)^ and allowing for /i oo. This expectation 
is supported by recent rigorous results [53] demonstrating that the dynamical equations 
known for the spherical spin glass model are faithfully reproduced in this limiting pro- 
cedure. Let us however note that the above choice goes beyond the immediate scope of 
the present paper as we decided to restrict our attention by considering only monoton- 
ically increasing confining potentials. Actually, investigating other classes is not at all 
a problem, and is certainly worth doing in view of the mentioned correspondence with 
the spherical model, where quite a few results on complexities of extrema of free energy 
landscapes were already reported in the literature, see for example a good review [51] and 
references therein. We relegate this issue, as well as a few other possible generalisations 
to subsequent publications. 
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APPENDICES 



A On a geometric origin of the identity (llUl) 



Let us start with recalling the well-known Poincare-Hopf index theorem. For any vector 
field F over any compact manifold M with no zeroes on the boundary dM holds the 
identity Ind{F) + Ind{dJF) = x{M), where xi^) is a topological invariant, the Euler 
characteristic of the manifold. The index Ind of a vector field F is the sum of indices 
of all the singular points, i.e. sum of indices signdet (i9jFj)|x=xfc corresponding to all 
isolated zeroes of the field: F(xfc) = . The field 9_F is a vector field defined in terms 
of the values of F on a subset of the boundary dM, such that the original field F points 
inward on this subset of the boundary. Precise definition of the boundary component of 
the formula is immaterial for our purposes. 

Consider a random surface H(x.) = ^x^A'x + V(x), with x G M^, V^(x) being a mean- 
zero random Gaussian function with the covariance Eq.Q, and the first term being a 
quadratic form involving an arbitrary non-singular real symmetric matrix A, with det A ^ 
0. Take any compact manifold M in such that (i) all the stationary points x^ of 
H[x) are with probability one belong to the interior of that manifold and (ii) along the 
boundary dM of the manifold the infiuence of random potential is already negligible (the 
simplest choice of M would be a ball |x| < L of a very large radius L —>■ oo). Then 
topological properties of the gradient field F = (9if(x)/(9x along the boundary dM are 
determined by the gradient of the first term, i.e. by the vector field F^ = Ax, so that 
Ind{d-F) = Ind{dJF a)- Then the Poincare-Hopf theorem implies Ind{dH{x.)/dx.) -t- 
Ind{d-FA.) = x{M). On the other hand, the field F^ has the only zero at the origin, 
so that Ind{FA) = sign det A, and the Poincare-Hopf theorem appied to Fa requires 
sign det A + Ind{d-FA) = x(M). Comparing these two relations we see that 

Ind{dH{^)/d^) = J2 sign det ( —-^ j U=x, = sign det A (99) 
k=i V ^ ^/ 

for every realization of the random surface H{x.). 

On the other hand, the Dirac's 5-functional measure satisfies the fundamental identity 

fdH\ , , / d^H^ 



which in fact underlies the Kac-Rice formula Eq(j4]). As above, the summation over x^ 
goes over all isolated zeros of the gradient dH/d^, and 5 (^) = Hili ^ (l^)' '^^^^ ^^^^ 
allows us to convert the sum in the left hand-side of (!99|) to an integral and rewrite (199!) 
as the identity 
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valid for any realization of the random potential V^(x). We therefore can average (llOip over 
the realizations, and take into account (i) independence of the first and second derivatives 
for Gaussian-distributed random functions, and (ii) stationarity of the random potential, 
Eq-dSD- In this way we come to the relation 

sign det A ^ (det (i + ^) ) ^ ft ) ) 

Using the standard Fourier integrals representation for the 5— functional factors, and 
performing the averaging over the gaussian gradients with the covariances ( cf. fl^T]) ) 

§■§-) = we see that 

'fr(i)) - /e---*-*'<»'l^ (103) 

iV „T ;t 

.g ^ (104) 



(27r|/'(0)|/iV)^A 

Substituting the last expression to Eq. fll02p and performing the Gaussian integral over 

r . ^-1-1/2 

yields the factor det A =1/1 '^^t A\. Then eq.f ll02p is reduced to the relation 

equivalent to the identity Eq. flTU]) . which is thus verified for non-singular real symmetric 
matrices A. By analytic continuation it is extended to singular well. 



B Distribution of the diagonal element of the GOE 
resolvent 

Our goal is to calculate the probability distribution V{G) of the diagonal element of the 
resolvent for GOE matrices in the large-N limit, for a given real value of s. 
Following the standard route we first evaluate the characteristic function 

X{p) = (exp{ipG^^(x)})GOE (106) 

Our first observation is that the result can depend on x only via the modulus |x| due to 
the rotational invariance of the GOE probability density. This implies that we can choose 
X = |x|e, with e = (1, 0, . . . , 0), so that 

H + si ~[ s + An 
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where A„, en stand for the eigenvalues and the corresponding eigenvectors of the GOE 
matrix H, and (e, e„) is the scalar product. Remembering that the eigenvectors en and 
the eigenvalues A„ of the GOE matrices are statistically independent, we perform the 
averaging over the eigenvectors first. To this end we recall that GOE eigenvectors are 
(i) mutually orthogonal and (ii) uniformly distributed over the unit sphere (e„, e^) = 1. As 
is well-known, these conditions imply that in the large- limit the projections Vn = (e, e„) 
behave like independent Gaussian variables with zero mean and variance {v'^)goe = 1/N. 
Denoting (. . .)„ the averaging over these variables, we can write for the characteristic 
function 



N 



{ex.p{ipGH{:>^)})v = JJ 



n=l 



N 

n 

n=l 



N 



S + Xn 



1/2 



dVr, 



N 



, ■exp{-^—vl + ip^'^ V} (108) 



N 



Remembering the remaining averaging over the JPD of GOE eigenvalues A„, we see that 
the characteristic function fll06p in the large-A^ limit can be written as the expectation 
value of the ratio of square roots of the characteristic polynomials 



x(p)U^oo = lim 

N^oo 



det(s + H) 



1/2 



det(s-2i^x2 + if) 



1/2 



(109) 



GOE 



Precisely that expectation value was already calculated earlier in a different context, see 
Eq.(35) and Eq.(48) of the paper ^7\, with the result 



x{p)\ 



exp 



-sp- \p\7njsc{s) 
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\s\ < 2 



;iio) 



The distribution V{G) immediately follows from this expression after the Fourier-transform. 
Defining G = (j'i^(x)/x^ we see that the quantity is Cauchy-distributed: 



V{G) 



7r2^|,(s) + (G-f) 



5iG-l 



\s\ < 2 
Isl > 2 



:iiii 



C An overview of the Dean-Majumdar functional in- 
tegral approach 

Our aim in this section is the calculation of large-A^ asymptotics of the required GOE 
averages. The most economic way of arriving to the desirable expressions known to us 
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relies on a heuristic method introduced by Dean and Majumdar ^6]. We however have 
every reason to beheve that one can arrive to the same results by employing a rigorous 
(and, necessarily, tedious) mathematical procedures described for a very closely related 
problem in a paper by A. Boutet de Monvel, L.Pastur and M. Scherbina 

The quantities of interest for us are the GOE averages featuring in equations fl35l) and 
dMD, i.e.: 



\det{H + si)\ 



GOE 



OC 



■ N 

n 

i=l 



N 



l<i<j<Af 



;ii2) 



i=l 



9{H + si) det{H + si) 



GOE 



OC 



■ N 

n 

.1=1 



dXi 



N 



l<i<j<N 



;ii3) 



i=l 



where 6'-factor is equal to unity when its argument is a positive definite matrix and is 
zero otherwise, and Aj stands for real eigenvalues of the matrix H. We have used that 
in each case the functions to be averaged depend only on the matrix eigenvalues. The 
average is taken over an ensemble of matrices with the probability density V{H)dH oc 
exp(— A^Trif^/4) invariant with respect to orthogonal transformations H —>■ OHO^, and 
this allows us to perform in a standard way [49] the integration over N{N — l)/2 angular 
variables (eigenvectors of H). The procedure yields just an overall normalization factor, 
and the remaining expression is given by flll2p . At this point it is worth mentioning, 
that the GOE averages of the type as above can be viewed as performed over a Gibbs 
measure describing an interacting gas of eigenvalues, and there are methods to developed 
a rigorous mean- field description of such a model in the large- A^ limit |55] . 

In a similar in spirit, but much less formal way Dean and Majumdar suggested to 
calculate integrals of this type by replacing the multiple integration over the eigenvalues 
Aj with a functional integration over the density of eigenvalues, defined as 



1=1 

Exploiting this density, we have the formal identities 

N 

J]0(A,) = expAT / ln0(A)pjv(A)rfA, 
1[^{K,X,) = exp— I lmlj{X,X')pN{X)pN{X')dXdX' 



;il4) 



(115) 
(116) 
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for a suitable choice of functions 0(A) and ■0(A, A') = '?/'(A', A). 

In this way the integrands of the multivariable integrals flll2p can be viewed as func- 
tionals /[p] of the density Pat (A), and introducing the (functional) Dirac 5— function, one 
can formally write for such type of integrals 



■ N 

n 

i=l 



N\ 



mm 



■ N 

n 

1=1 



5[p- Pn] 



;il7) 



where /[p] is the functional / expressed in terms of the density p. The integrals over the 
eigenvalues on the right hand side of this expression can be conveniently evaluated after 
using the standard formal Fourier representation for the functional 5— function: 



^[p- Pn] 



V[uj] exp i / uj{t)[p{t) - pN{t)]dt 



f roc ^ \ 

T> [uj] exp \ iN uj {t)p{t)dt - i^uj{\i) 



;il8) 

;ii9) 



with a suitably normalized measure D[ct;]. The integration over the eigenvalues is now 
trivially performed, resulting in 



n 

i=l 



dXi 



5[p - Ptv] oc J V[u;'] exp (^iN j uj\t)p{t)dt + Nln(^j^ 



(120) 

For ^ 1 the main contribution to the above integral should come from the stationary 
point of the exponent with respect to variations in the field uj' satisfying the equation 



P(A) 



exp(— ^(^^'{A)) 
/>xp(-i<^'(A'))dA' 





if a < A < 6, 
otherwise. 



(121) 



It is clear from the above that p{X)dX = 1, and that —iuj'{X) = Inp(A) + const for 
a < X < b. This yields, up to an overall constant factor, the relation 



■ N 

n 

i=l 



dXi 



5[p - Pn] oc exp (^-N j p{t) \n{p{t))dt^ 



;i22) 



Note, that looking at our system as a kind of "eigenvalue gas", this factor is simply a 
standard entropic contribution associated with the density of particles. 

Let us apply this method for our problem, i.e. to the calculation of the required GOE 
averages. Denoting the range of integration in each case by TZg and TZm-, we have for 
i e {s,m} 



j P[p]exp (^-^g^[p]+Nr,[p] 



(123) 
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where 

GAn\ = ^ 



SAp] = i / AV(A)rfA- / / /5(A)p(A')ln(|A-A'|)rfAc/A', (124) 



Tj[p\ = [ p(A)ln(|A + s|)rfA- [ p(A) lnp(A)dA. (125) 

In the hmit iV — >■ cxd in (11231) the main contribution to the functional integral in (11231) 
comes obviously from the value of p which minimises the functional Qj. The stationary 
condition is found in the standard variational procedure after incorporating the normali- 
sation condition p{X)dX = 1 via a Lagrange multiplier. The resulting integral equation 
reads: 

\2 r 

+ C= p(A')ln(|A- A'|)dA' for A G 7^J• . (126) 
4 Jn, 

Differentiating this equation with respect to A gives 

A f piX')dX' , , ^ 

where the integral must now be understood as a Cauchy principal value. 

Solution of singular integral equations of the type ( 1127^ is discussed extensively in [58] . 
The necessary inversion formula essentially depends on whether the solution is required to 
be bounded at each end. In general, given a function g{X) for A G (a, h) and the integral 
equation 

-' f{\')d\' 



■ A -A' 

there is a unique solution /(A) which remains bounded at the endpoints a and h provided 
that holds 

g{X')dX' Q _ 



la ^{b-X'){X' -a) 
In our case, the above condition implies 

= f , + (128) 

Hence, such a solution only exists in the case a = —b = L/2, and in the latter case is 
given by the inversion formula 



JLVA - A2 /-^/^ X'dX' , , 

p A = ^ — '— / , , 129 

^ 27r2 v^LV4-A'2(A' - A) ' ^ ^ 



Vi^V4 - A2 
27r2 

V^V4 - A2 
27r 



dX! ^ /-^/^ ^y, 

+ A 



L/2 ^L2/4 - A'2 y_L/2 ^L2/4- A'2(A' - A) 



, (130) 
(131) 
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where in the final fine we have used the identity 

/ w w ; = for.G(a,6). (132) 

J a \/[o — X) [x — a)[x — Zj 

The normalization condition fixes L, and the resulting eigenvalue density is given by 

PscW = T^V4- A2, for A G [-2, 2] (133) 

ZTT 

and zero otherwise, which is just the well-known Wigner semi-circle law. In such a case 
the only dependence on the variable s in the exponent of 'Ds{s) comes from the first term 
/^2Psc(A) ln(|A + s\)dX in fll25p which leads to Eq.(37) as described in the text. 

At the same time, when evaluating Vm{s) we have an additional constraint on the 
density p, as the latter must vanish for A < —s. Obviously, the solution fll33p can satisfy 
such a constraint only as long as s > 2, and has to be modified in the opposite case s > 2. 
As shown in [SB] there always exists a solution of this type of integral equations which is 
bounded only at one end of the integration range. The solution which remains bounded 
at the upper end of the integration range is given by 



Pdm(A) = —A , , ^ / \\ ^ -, (134) 



27r2 V A + s V L - s - A' A - A 

where L is a constant to be determined by the normalization. The integral was further 
evaluated by Dean and Majumdar and the resulting density is given, for our choice of the 
GOE measure, by 




L-\-s 



[L + 2A] if < A + s < L, 



where this time 



PDA/ (A) = V ' " ^ - - . ^ - ^-^35) 

] otherwise, 

L = ^ + Vs2 + 12) . (136) 

Note that s — 2 implies L ^ 4 and the Dean-Majumdar density ( 1135^ reverts in this 
limit to the Wigner semicircular law f ll33p . 

Inserting the above function pdm in the definitions of Qm and 7^, see f ll24p .( 11251) and 
performing the integrations with help of Mathematica yields 

Qm\pDwi\ = G^,^{s) = (72s^ - - 30s V12 + s2 - s V12 + s2 + 54(3 + 4 In 6) 

210 v 

- ln(s + Vs^ + 12) , (137) 

which is our expression Eq.(69), as well as a formula for %n[pDAi]- The latter formula is 
however irrelevant for finding the complexity of minima, as our analysis reveals that the 
complexity is determined by value s 2 when obviously %n[pDAi] %n[Psc]- 
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D Analysis of the replica-symmetric solution for a 
spherically-symmetric confining potential. 



The calculation of the de-Almeida- Thouless condition in the general energy surface of 
the form ffTTl) requires only very minor modifications in comparison with the case of a 
parabohc confinement ([6]) considered in detail in [13] . Applying the procedure of Ref . [13] 
yields the following exact expression for the averaged replicated partition function: 



{ZD = C^,„Ar^"/2e'^^"^W / (detQ)-("+^)/' e-^^*"(«) dQ 

Jq>q 



:i38) 



where 



n ^ 



-{qaa + qbb) - Qab 



(139) 



a=l a<b 

where Cjy^n is a known numerical constant and N is assumed to satisfy the constraint 
N > n. 

The form of the integrand in Eq. fll38p is precisely one required for the possibility of 
evaluating the replicated partition function in the limit — oo by the Laplace ("saddle- 
point") method. The free energy is then given by 



hm Lf) = -|ln(27re) - ^/(O) + lim i$„(Q) 



(140) 



where the entries of the matrix Q are chosen to satisfy the stationarity conditions: 
^ ^'^^ = for a < b. This yields, in general, the system of n{n + l)/2 equations: 



1 



b{^a) 



Maa + Qbb) - qab 



0, a = 1,2, 



n 



(141) 



and 



--p[Q-\b + Pf' 



^{qaa + qbb) - qab 



0, a^b 



(142) 



where f'{x) stands for the derivative df /dx and we used the convention ( |T2l) . 

The Replica Symmetric Ansatz amounts to searching for a solution to Eqs. (1141111421) 
within subspace of matrices Q = Qrs > such that qaa = qd, for any a = 1, . . .n, and 
qa<b = qo- The system of equations is easy to solve and to obtain in the replica limit 
n — i> the following relations: 



qd = 2/ — 



qo 



(143) 



where we denoted for brevity fid = fi{qd 
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The stability analysis of this solution amounts to expanding the function $„(Q) in 
Eq. fll39p around the extremum point up Q = Qrs to the second order in deviations: 
$ = (^gp -)- (5$ -)- The stationarity condition amounts to 5$ = yielding the 

system f ll41l) - fll42p . The term is a quadratic form in independent fluctuation variables 
5qab,a < b and can be generally written as 5^$ = X](a6),{cd) ^(l{ab)G(ab),{cd) Sq^cd)- As usual 
the stable extremum corresponds to the positive definite quadratic form, and along the 
critical line the quadratic form becomes semi-definite. Checking positive definiteness of 
(5^<I> amounts to finding the (generalized) eigenvalues A of the matrix G(^ab),{cd) 
It is easy to see that in our case the quadratic form can be written as 

S'^ = ^ j^^habf + ^Tr [5Q{QRsr'mQRsr'] -^Y. f"(Dab)SDl , (144) 

a ab 

where we introduced short-hand notations 5Dab = \{5qaa + ^qbb — '^^qab) and Dab = qd — qab- 
The corresponding eigen-equations for generalized n{n + l)/2 component eigenvectors rji^ab) 
with a < b can be written straightforwardly by repeating the analysis of [I3] and are given 
for T > by: 



dqabdqcd ' 



E iQ~Rs)ac^ic^) {QRs)db+Y^f"(Dab){5Dab)+S^, 

cd 



;^/i'('?d)-$^r(/^ae)(5I^ac) 



A^ 



V{ah) 



(145) 

As the entries of the matrix Qrs are given by qab = 90+ (o'd " %)^ab) its inverse Q ^ has 
the same form {Q~^)ab = Po + {Vd—Pol^abi with po and Pd—Po are given in the limit n 
by 

= —/ ^ = ^/ ~ ' Pd-Po = = 7^ (146) 

[(id-q^Y \fidj qd-qo T 

Now we can follow faithfully the lines of the classical work by De Almeida and Thouless 
[3] and to provide an explicit construction of the families of eigenvectors with components 
T](^ab) of different symmetry. For the model with parabolic confinement this construction 
was discussed in [13j and goes through here without any modification. There are three 
different families of eigenvectors, first two yielding only eigenvalues A* with positive real 
part, hence stable. The dangerous third family of eigenvectors is that satisfying the 
constraints: 

d 

The equations Eq. fll45p are then reduced to a single equation 



(147) 



iPd-Po)^ - T^figd-qo) 



%b) 



(148) 



and substituting here (11461) we find that the replica symmetric solution is stable as long 
as 

/ T 



(149) 
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with the equahty fx^ = f" providing the condition of the rephca symmetry breaking 

transition. Solving the latter together with Eg. (11431) is easy after introducing the auxilliary 
variable ^ = This gives finally the transition temperature line in the form 

Tat = ry^rF), (150) 

where r satisfies the equation 

/^(^-;^) = v7^- (151) 

In particular, the transition at zero temperature T ^ requires t —>■ 0, hence the 
replica symmetric solution at zero temperature is stable as long as the inequality 

/^(-||y)>v^- (152) 

is satisfied, with the equality sign standing for the transition condition to the region with 
broken replica symmetry. 



E Complexity of stationary points with a given index 

In this appendix we outline the calculation of complexity for the stationary points with a 
given index T{H) (the number of negative eigenvalues of the Hessian H) for our model. 
We closely follow the method developed by Bray and Dean in [51], which allows one to 
perform the calculation for the extensive values of the index scaled with N a.s I = aN, 
where < a < 1. For simplicity we consider only the case of parabolic confinement, with 
constant fi{z) = /i. 

Repeating the same steps as in ^J, the expected value of the number Af{fi,a) of 
stationary points with a given value of a can be straightforwardly shown to be given by 
(cf. Eq.©): 

s^2^y N + 2 V TT j 



X 



J dH\ <\ei{slN + H)\5 {AsIn + H) - No) g-f [^^'-iv^(Tri?)^ 



where we denoted s = /i/ a/ /"(O) = fJ-/ fJ-cr- 

The integration over the real, symmetric matrix H is in the standard way reduced to 
an integral over its eigenvalues A, , and further to the functional integral over the mean 
eigenvalue density, Eq. (lll4p . see Appendix B. This gives: 

(A/'(/^,«)) = ^%^, (154) 
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where 



X(5(^y" p{\)d\ - 5 (^j p{\)d\-a]. (155) 
The expressions appearing in the exponential in (11551) are given by 



oo noo 



p{\)p{\')\Yi{\\-\'\)d\d\', (156) 



oo J — oo 



and 



/oo /»oo 
p(A)ln(|s + A|)rfA- / p(A)hip(A)c/A 
oo J —oo 

The denominator D appearing in (I154p is a normahzation factor which can also be rep- 
resented functional integral and is given by 

D = j V[p] exp{~N^S2[p] + NS[[p])6 (^j p{\)d\ - 1^ , (158) 

where ^ 

S[[p] = - 1^ p{X) In p{X)dX~^^^^^(^l^ p{X)XdX^ . (159) 

As — s> oo, the main contribution to the functional integral in (11551) comes from the 
function p which minimizes iS2[p]. By using the definition of the mean eigenvalue, 



oo 



A := / p(A)ArfA, (160) 



oo 



and further introducing the function / defined as 

/(a;)=p(x + A), (161) 
we rewrite S2[p] in the following form: 

-| POO -1 POO poo 

S2[f] = j fix)x'dx-- / fix)fix')\ni\x-x'\)dxdx'. (162) 

^ J —oo ^ J —oo J —oo 
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We now need to find the function / wliicli minimizes (11621) subject to tlie constraint 
J f{x)dx = 1. A simple variational calculation (cf. Appendix C) sliows tliat tlie minimizer 
is tlie usual Wigner semi-circle density, given by 



/(^) = (163) 
Zn 

The value of A can now be fixed using the restriction on the index given by the second 
(5— functional factor in (I155p . This leads to 

2 /"^^/^ 



TT 



a = - Vl - x^dx, (164) 



1 



where h = s + \. Each value of a G [0, 1] then corresponds to a unique value of h E [0, 2]. 

We can apply the same variational procedure to the normalization factor D. The 
only difference is that the functional integration in the denominator does not contain 
the restriction on index, which results in second term in S[ minimized by the value 
A = 0. On the other hand, the terms containing ^2 appear in both D and and 
therefore cancel from their ratio in (11541) when we apply the asymptotic evaluation of the 
integrals in the limit —>■ oo. Likewise, the term of the form J p(A) log p{X)dX will cancel 
between the numerator and the denominator. Consequently, the relation (11541) is reduced 
asymptotically to 



{Af{fi,a)) ~ l^exp(Ar 



' f{x)H\x + h\)-h' 

-2 



(165) 



where A = h — s, and h is related to a by Eq. (11641) . For the semicircular form of f{x) the 
integral above can be explicitly calculated: 

1 /•2 ^ 1 

— V4-x2 \n\x + h\dx = — + — . (166) 
2nJ_2 2 4 

Substituting this result to Eq. (I165P yields the final form of the complexity corresponding 
to a given index: 

S(.,a) = -l + ^-^^^-ln.. (167) 

The analysis of this expression is convenient to perform separately for a = and 
a > 0. 

1. The relation (11641) for a = implies h = 2. Inserting this in (I167P gives 

1 (2- s)"^ 

ns,0) = --^-^-\ogs. (168) 
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This expression is zero at s = 1, and taking the derivative with respect to s gives 



^S(s,0) =2- + <0 (169) 

with equahty achieved only for s = 1. We see that S(s, 0) is decreasing with s and 
hence is positive for < s < 1. Also ^S(s, 0) = ^ - 1, ^S(s, 0) = Thus, 
the first non-vanishing derivative at s = 1 is the third derivative, implying that the 
complexity of stationary points with a = must vanish cubically as s — > 1. This 
agrees with the general analysis of complexity of minima performed earlier in the 
paper. 

2. Eq. fll64p implies h < 2 for a > 0. Differentiating (I167P with respect to s gives 

^S(s,a) = /i- ( s + - ) < 0. (170) 
OS \ s J 

Thus, S(s,a) is strictly decreasing as a function of s, and vanish linearly at some 
critical point Scr which can be easily found numerically as a function of a. Observing 
that 

E(l,a) = -i^^<0. (171) 

we conclude that the critical value Scr(a) must satisfy Scr(«) < 1 for any a > 0. In 
other words, the complexity of stationary points with any extensive index I = 0{N) 
is still negative at the point of ergodicity breaking s = 1 (i.e. /i = ficr = a/ /"(O) 
in the original notations), and starts to be positive already inside the phase with 
broken ergodicity: /icr(« > 0) < a/ /"(O). 

If we considered the annealed average of the total number of stationary points, rather 
than those with a fixed index, this would be equivalent to integrating (A/'(/i, a)) over all 
values of a. In the limit ^ 1 such an integral will obviously be dominated by the value 
of a which maximizes E(s, a). As da/dh ^ for \h\ < 2 the maximum occurs when 

= |^S(.,a) = .-^, (172) 
so that h = 2s. At this point of maximum 

^tot{s) = ^Y^-\ogs. (173) 

This shows that complexity of all stationary points tends to zero as s ^ 1. and by taking 
derivatives we find that at s = 1 the first derivative vanishes but the second one is non- 
zero. Consequently, the (annealed) complexity related to the total number of stationary 
points vanishes as (1 — s)^ as s — 1, as indeed was found in |18] by a different method. 
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Moreover, the last result can be used to show that the critical value Scr{ct) , - defined 
as s at which the complexity with a given a vanishes, - must monotonously decrease with 
a increasing. First observe that in fact we have shown that E(/i/2,a) > which implies 
Scr(a) > h/2. The value Scrih) by definition solves S(scr(/i), h) = 0. As /i decreases with 
a increasing, it is enough to show that dscr/dh > 0. Taking the derivative with respect 
to h and rearranging gives 

dscr ^ _( dE/dh\ ^ Scr -h/2 

dh \dJ:/ds J s,r + S-^-h ■ ^ ' 

as required. 
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